Визуализация биомедицинских данных. Д/з №2
Anastasia Soroka
2023-11-12
1. Импорт данных и первичный анализ
data <- readRDS("life_expectancy_data.RDS")Описание данных
data_info <- describe(data, na.rm = TRUE, skew = FALSE, ranges = TRUE)
data_info## vars n mean sd
## Country* 1 195 9.800000e+01 5.644000e+01
## Year 2 195 2.019000e+03 0.000000e+00
## Gender* 3 195 1.000000e+00 0.000000e+00
## Life expectancy 4 195 7.552000e+01 7.690000e+00
## Unemployment 5 195 8.600000e+00 6.970000e+00
## Infant Mortality 6 195 1.961000e+01 1.693000e+01
## GDP 7 195 4.659960e+11 1.933272e+12
## GNI 8 195 4.864381e+11 1.956716e+12
## Clean fuels and cooking technologies 9 195 6.598000e+01 3.632000e+01
## Per Capita 10 195 1.682098e+04 2.444060e+04
## Mortality caused by road traffic injury 11 195 1.706000e+01 1.037000e+01
## Tuberculosis Incidence 12 195 1.037500e+02 1.342800e+02
## DPT Immunization 13 195 8.799000e+01 1.240000e+01
## HepB3 Immunization 14 195 8.676000e+01 1.273000e+01
## Measles Immunization 15 195 8.731000e+01 1.319000e+01
## Hospital beds 16 195 3.000000e+00 2.370000e+00
## Basic sanitation services 17 195 7.738000e+01 2.777000e+01
## Tuberculosis treatment 18 195 7.757000e+01 1.705000e+01
## Urban population 19 195 5.912000e+01 2.329000e+01
## Rural population 20 195 4.088000e+01 2.329000e+01
## Non-communicable Mortality 21 195 1.705000e+01 7.080000e+00
## Sucide Rate 22 195 4.800000e+00 3.780000e+00
## continent* 23 195 2.670000e+00 1.310000e+00
## min max range
## Country* 1.00 1.950000e+02 1.940000e+02
## Year 2019.00 2.019000e+03 0.000000e+00
## Gender* 1.00 1.000000e+00 0.000000e+00
## Life expectancy 55.49 8.810000e+01 3.260000e+01
## Unemployment 0.18 3.644000e+01 3.626000e+01
## Infant Mortality 1.40 7.580000e+01 7.440000e+01
## GDP 188391770.64 2.143322e+13 2.143304e+13
## GNI 375392118.22 2.170865e+13 2.170827e+13
## Clean fuels and cooking technologies 0.00 1.000000e+02 1.000000e+02
## Per Capita 228.21 1.758139e+05 1.755857e+05
## Mortality caused by road traffic injury 0.00 6.460000e+01 6.460000e+01
## Tuberculosis Incidence 0.00 6.540000e+02 6.540000e+02
## DPT Immunization 35.00 9.900000e+01 6.400000e+01
## HepB3 Immunization 35.00 9.900000e+01 6.400000e+01
## Measles Immunization 37.00 9.900000e+01 6.200000e+01
## Hospital beds 0.20 1.371000e+01 1.351000e+01
## Basic sanitation services 8.63 1.000000e+02 9.137000e+01
## Tuberculosis treatment 0.00 1.000000e+02 1.000000e+02
## Urban population 13.25 1.000000e+02 8.675000e+01
## Rural population 0.00 8.675000e+01 8.675000e+01
## Non-communicable Mortality 4.40 4.370000e+01 3.930000e+01
## Sucide Rate 0.30 3.010000e+01 2.980000e+01
## continent* 1.00 5.000000e+00 4.000000e+00
## se
## Country* 4.040000e+00
## Year 0.000000e+00
## Gender* 0.000000e+00
## Life expectancy 5.500000e-01
## Unemployment 5.000000e-01
## Infant Mortality 1.210000e+00
## GDP 1.384445e+11
## GNI 1.401234e+11
## Clean fuels and cooking technologies 2.600000e+00
## Per Capita 1.750230e+03
## Mortality caused by road traffic injury 7.400000e-01
## Tuberculosis Incidence 9.620000e+00
## DPT Immunization 8.900000e-01
## HepB3 Immunization 9.100000e-01
## Measles Immunization 9.400000e-01
## Hospital beds 1.700000e-01
## Basic sanitation services 1.990000e+00
## Tuberculosis treatment 1.220000e+00
## Urban population 1.670000e+00
## Rural population 1.670000e+00
## Non-communicable Mortality 5.100000e-01
## Sucide Rate 2.700000e-01
## continent* 9.000000e-02na_count_cols <- colSums(is.na(data))
na_count_cols #Сколько N/A в каждом столбце## Country Year
## 0 0
## Gender Life expectancy
## 0 0
## Unemployment Infant Mortality
## 0 0
## GDP GNI
## 0 0
## Clean fuels and cooking technologies Per Capita
## 0 0
## Mortality caused by road traffic injury Tuberculosis Incidence
## 0 0
## DPT Immunization HepB3 Immunization
## 0 0
## Measles Immunization Hospital beds
## 0 0
## Basic sanitation services Tuberculosis treatment
## 0 0
## Urban population Rural population
## 0 0
## Non-communicable Mortality Sucide Rate
## 0 0
## continent
## 0
2. Интерактивный plotly график для двух нумерических переменных (), раскрашенный по колонке континента, на котором расположена страна
plot_ly(
data[data$`Life expectancy` != 0,],
x = ~ `Life expectancy`,
color = ~continent,
type = "box"
)plot_ly(
data[data$`Sucide Rate` != 0,],
x = ~ `Sucide Rate`,
color = ~continent,
type = "box"
)
3. Cравнение распределений колонки Life expectancy между группами стран Африки и Америки. Вид статистического теста определите самостоятельно. Визуализируйте результат через библиотеку rstatix.
data_filtered <- data %>%
filter(continent %in% c("Africa", "Americas"))
test_data_filtered <- data_filtered %>%
t_test(`Life expectancy` ~ continent,
data = .,
paired = FALSE,
var.equal = TRUE)
print(test_data_filtered)## # A tibble: 1 × 8
## .y. group1 group2 n1 n2 statistic df p
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl>
## 1 Life expectancy Africa Americas 52 38 -11.4 88 4.12e-19Поскольку p-value ниже заданного значения 0.05, то мы отвергаем Н0 и можем утверждать, что продолжительность жизни отличается в странах Африки и Америки.
sum_stats <- data_filtered %>%
get_summary_stats(`Life expectancy`, type = "mean_sd")
sum_stats## # A tibble: 1 × 4
## variable n mean sd
## <fct> <dbl> <dbl> <dbl>
## 1 Life expectancy 90 71.4 8.05library(ggpubr)
ggqqplot(data_filtered[data_filtered$`Life expectancy`, ],
x = "Life expectancy", facet.by = "continent")
#Данные не выглядят нормально распределенными, используем тест Манна-Уитни вместо Т-тестаstat.test <- data_filtered[data_filtered$`Life expectancy`, ] %>%
wilcox_test(`Life expectancy` ~ continent) %>%
add_xy_position(x = "Life expectancy")
stat.test## # A tibble: 1 × 11
## .y. group1 group2 n1 n2 statistic p y.position groups xmin
## <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <name> <dbl>
## 1 Life exp… Africa Ameri… 47 43 72 3.25e-14 84.1 <chr> NA
## # ℹ 1 more variable: xmax <dbl>library(rstatix)
ggboxplot(
data_filtered[data_filtered$`Life expectancy`, ],
x = "continent", y = "Life expectancy",
ylab = "Life expectancy", xlab = "Continent",
add = "jitter"
) +
labs(subtitle = get_test_label(stat.test, detailed = FALSE)) +
stat_pvalue_manual(stat.test, tip.length = 0) 
4. Сделайте новый датафрейм, в котором оставите все численные колонки кроме Year. Сделайте корреляционный анализ этих данных. Постройте два любых типа графиков для визуализации корреляций.
data_filtered_2 <- data %>%
select_if(is.numeric) %>%
select(-Year, -`Life expectancy`)
data_f2_cor <- cor(data_filtered_2) library(corrplot)## corrplot 0.92 loadedcorrplot(data_f2_cor, method = 'number', number.cex = 0.8, tl.cex = 0.7)
corrplot(data_f2_cor, method = "color", type = "lower",
addCoef.col = "grey29", diag = FALSE,
cl.pos = "r", tl.col = "grey11",
col = COL2('RdBu', 10),
tl.cex = 0.8)
5. Постройте иерархическую кластеризацию на этом датафрейме.
library(factoextra)## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa#Стандартизация данных
data_scaled <- scale(data_filtered_2)
#Матрица дистанций
distances_data_scaled <- dist(data_scaled,
method = "euclidean"
)
as.matrix(distances_data_scaled)[1:6, 1:6]## 1 2 3 4 5 6
## 1 0.000000 7.390287 6.144642 4.404280 6.468405 7.723136
## 2 7.390287 0.000000 2.610829 7.637350 3.346255 3.630919
## 3 6.144642 2.610829 0.000000 6.049814 4.350330 3.456356
## 4 4.404280 7.637350 6.049814 0.000000 7.886278 6.853670
## 5 6.468405 3.346255 4.350330 7.886278 0.000000 4.960145
## 6 7.723136 3.630919 3.456356 6.853670 4.960145 0.000000#Дендрограмма кластеров
data_scaled_hc <- hclust(d = distances_data_scaled,
method = "ward.D2")
#Визуализация
fviz_dend(data_scaled_hc,
cex = 0.1) # cex() - размер лейблов
6. Сделайте одновременный график heatmap и иерархической кластеризации. Содержательно интерпретируйте результат.
library(pheatmap)
plot3 <- pheatmap(data_scaled,
show_rownames = FALSE,
clustering_distance_rows = distances_data_scaled,
clustering_method = "ward.D2",
cutree_rows = 5,
cutree_cols = length(colnames(data_scaled)),
angle_col = 45,
main = "Heatmap + hierarchical clustering")
plot3
В данном графике дендрограмма сверху показывает иерархическую кластеризацию переменных. Кластеры отображают какие переменные схожи, а высота соединения показывает степень различия между группами переменных. Два больших кластера. Первый с переменными, связанными с уровнем безработицы, туберкулезом и его лечением, разными Mortality. Во втором кластере переменные GDP и GNI являются похожими и стоят отдельно от всех остальных переменных в этом кластере. Переменные связанные с иммунизацией сгруппировались в один схожий кластер. Вертикальная дендрограмма (сбоку) отражает иерархическую кластеризацию наблюдений (объектов). Здесь выделяется 5 групп (одна очень маленькая с высокими стандартизированными значениями переменных GDP, GNI). Эти 5 групп, вероятно, более-менее отражают группы по континентам.
counts <- table(data$continent)
# Вывод результатов
print(counts)##
## Africa Americas Asia Europe Oceania
## 52 38 42 48 15
7. Проведите PCA анализ на этих данных. Проинтерпретируйте результат.
library(FactoMineR)
data_full.pca <- prcomp(data_filtered_2,
scale = T) summary(data_full.pca)## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 2.6027 1.4816 1.3933 1.16723 1.08021 0.96332 0.92844
## Proportion of Variance 0.3764 0.1220 0.1078 0.07569 0.06483 0.05155 0.04789
## Cumulative Proportion 0.3764 0.4983 0.6061 0.68184 0.74666 0.79822 0.84611
## PC8 PC9 PC10 PC11 PC12 PC13 PC14
## Standard deviation 0.84336 0.69074 0.68425 0.58700 0.54757 0.42231 0.3550
## Proportion of Variance 0.03951 0.02651 0.02601 0.01914 0.01666 0.00991 0.0070
## Cumulative Proportion 0.88562 0.91213 0.93814 0.95728 0.97394 0.98385 0.9909
## PC15 PC16 PC17 PC18
## Standard deviation 0.34451 0.20224 0.07146 3.138e-16
## Proportion of Variance 0.00659 0.00227 0.00028 0.000e+00
## Cumulative Proportion 0.99744 0.99972 1.00000 1.000e+00fviz_eig(data_full.pca, addlabels = T, ylim = c(0, 40))
Первые две компоненты объясняют около 50% дисперсии.
fviz_pca_var(data_full.pca, col.var = "contrib", labelsize = 3)
fviz_pca_var(data_full.pca,
select.var = list(contrib = 5), # Задаём число здесь
col.var = "contrib")
fviz_contrib(data_full.pca, choice = "var", axes = 1, top = 24) # 1
fviz_contrib(data_full.pca, choice = "var", axes = 2, top = 24) # 2
fviz_contrib(data_full.pca, choice = "var", axes = 3, top = 24) # 3
Infant Mortality and Clean fuels and cooking technologies стали самыми важными переменными с точки зрениях их вариации в PC1. Иммунизации стали тремя самыми важными переменными с точки зрениях их вариации в PC2.
8.Постройте biplot график для PCA. Раскрасьте его по значениям континентов. Переведите его в plotly. Желательно, чтобы при наведении на точку, вы могли видеть название страны.
ggbiplot(data_full.pca,
scale=0, alpha = 0.1) +
theme_minimal()
Более осмысленным biplot становится при использовании кластерных методов, с помощью которых мы можем разделить наблюдения на группы. Посмотрим, наблюдается ли разница между группами по continents:
# Сделаем корректные данные для группировки по continents.
data_filtered_2_with_continent <- data_filtered_2 %>%
mutate(continent = data$continent)
# Визуализируем с группировкой по diabetes (для этого переменную нужно сделать фактором)
ggbiplot(data_full.pca,
scale=0,
groups = as.factor(data_filtered_2_with_continent$continent),
ellipse = T,
alpha = 0.2) +
theme_minimal()
pca_df <- as.data.frame(data_full.pca$x)
pca_df$country <- data$Country
pca_df$continent <- as.factor(data_filtered_2_with_continent$continent)
p <- plot_ly(pca_df, x = ~PC1, y = ~PC2, text = ~country, color = ~continent,
type = 'scatter', mode = 'markers',
marker = list(opacity = 0.7)) %>%
layout(title = 'PCA Plot',
xaxis = list(title = 'PC1'),
yaxis = list(title = 'PC2'))
p
9.Дайте содержательную интерпретацию PCA анализу.
Африка сильно отличается от Европы, Америки и Азии. Океания самая самая “разбросанная”, но в ней и меньше всего объектов. Infant Mortality and Clean fuels and cooking technologies стали самыми важными переменными с точки зрениях их вариации в PC1, а три вида иммунизации стали тремя самыми важными переменными с точки зрениях их вариации в PC2.
10. Сравните результаты отображения точек между алгоритмами PCA и UMAP.
Tidymodels approach
#install.packages("embed")
library(tidymodels)
library(embed)
umap_prep <- recipe(~., data = data_filtered_2) %>% # "техническая" строка, нужная для работы фреймворка tidymodels
step_normalize(all_predictors()) %>% # нормируем все колонки
step_umap(all_predictors()) %>% # проводим в UMAP. Используем стандартные настройки. Чтобы менять ключевой параметр (neighbors), нужно больше погружаться в машинное обучение
prep() %>% # "техническая" строка, нужная для работы фреймворка tidymodels. Мы выполняем все степы выше
juice() # Финальная строка - приводим результаты UMAP к стандартизированному датасетуВизуализиуем два первых измерения UMAP и добавим информацию о возрастных группах и диабет-статусе:
umap_prep %>%
ggplot(aes(UMAP1, UMAP2)) + # # можно добавить раскраску
geom_point(aes(color = as.character(data_filtered_2_with_continent$continent)),
alpha = 0.7, size = 2) +
labs(color = NULL) 
Отображение точек в целом схоже в PCA и в UMAP.
11. Давайте самостоятельно увидим, что снижение размерности – это группа методов, славящаяся своей неустойчивостью. Удалите 5 случайных колонок. Проведите PCA анализ. Повторите результат 3 раза. Наблюдаете ли вы изменения в куммулятивном проценте объяснённой вариации? В итоговом представлении данных на биплотах? С чем связаны изменения между тремя PCA?
set.seed(123)
columns_to_remove_1 <- sample(names(data_filtered_2), 5)
data_filtered_3 <- select(data_filtered_2, -all_of(columns_to_remove_1))
columns_to_remove_1## [1] "Urban population" "Tuberculosis treatment" "GDP"
## [4] "HepB3 Immunization" "Infant Mortality"columns_to_remove_2 <- sample(names(data_filtered_2), 5)
data_filtered_4 <- select(data_filtered_2, -all_of(columns_to_remove_2))
columns_to_remove_2## [1] "Measles Immunization"
## [2] "Clean fuels and cooking technologies"
## [3] "GNI"
## [4] "Tuberculosis treatment"
## [5] "Per Capita"columns_to_remove_3 <- sample(names(data_filtered_2), 5)
data_filtered_5 <- select(data_filtered_2, -all_of(columns_to_remove_3))
columns_to_remove_3## [1] "Clean fuels and cooking technologies"
## [2] "DPT Immunization"
## [3] "Basic sanitation services"
## [4] "GDP"
## [5] "Tuberculosis Incidence"data_full.pca_3 <- prcomp(data_filtered_3,
scale = T)
fviz_eig(data_full.pca_3, addlabels = T, ylim = c(0, 40))
data_full.pca_4 <- prcomp(data_filtered_4,
scale = T)
fviz_eig(data_full.pca_4, addlabels = T, ylim = c(0, 40))
data_full.pca_5 <- prcomp(data_filtered_5,
scale = T)
fviz_eig(data_full.pca_5, addlabels = T, ylim = c(0, 40))
Кумулятивный процент объясненной вариации для PC1 и PC2 при выбрасывании трех случайных переменных составил 50.8%, 54.2% и 47.6% (по сравнению с 49.8% при предыдущем анализе)
library(cowplot)
pl3 <- ggbiplot(data_full.pca_3, scale=0, alpha = 0.1) + theme_minimal()
pl4 <- ggbiplot(data_full.pca_4, scale=0, alpha = 0.1) + theme_minimal()
pl5 <- ggbiplot(data_full.pca_5, scale=0, alpha = 0.1) + theme_minimal()
plot_grid(pl3, pl4, pl5, ncol = 3)
При выбрасывании 5 переменных достаточно сильно изменилось отображение биплотов, хотя и некоторые переменные остались примерно в том же направлении (например, Hospital beds).
12. Давайте самостоятельно увидим, что снижение размерности – это группа методов, славящаяся своей неустойчивостью. Создайте две дамми-колонки о том: (1) принадлежит ли страна к африканскому континенту, (2) Океании. Проведите PCA вместе с ними, постройте биплоты. Проинтерпрейтируйте результат. Объясните, почему добавление дамми-колонок не совсем корректно в случае PCAнашего типа.
data_filtered_6 <- data_filtered_2_with_continent %>%
mutate(`is Africa` = as.integer(continent == "Africa"),
`is Oceanina` = as.integer(continent == "Oceania")) %>%
select(-continent)data_full.pca_6 <- prcomp(data_filtered_6,
scale = T)
fviz_eig(data_full.pca_6, addlabels = T, ylim = c(0, 40))
47.2% вариативности объясняют PC1 и PC2.
pl6 <- ggbiplot(data_full.pca_6, scale=0, alpha = 0.1) + theme_minimal()
pl6
# Сделаем корректные данные для группировки по continents.
data_filtered_6_with_continent <- data_filtered_6 %>%
mutate(continent = data$continent)
# Визуализируем с группировкой по diabetes (для этого переменную нужно сделать фактором)
ggbiplot(data_full.pca_6,
scale=0,
groups = as.factor(data_filtered_6_with_continent$continent),
ellipse = T,
alpha = 0.2) +
theme_minimal()
fviz_pca_var(data_full.pca_6, col.var = "contrib", labelsize = 3)
fviz_pca_var(data_full.pca_6,
select.var = list(contrib = 5), # Задаём число здесь
col.var = "contrib")
fviz_contrib(data_full.pca_6, choice = "var", axes = 1, top = 24) # 1
fviz_contrib(data_full.pca_6, choice = "var", axes = 2, top = 24) # 2
fviz_contrib(data_full.pca_6, choice = "var", axes = 3, top = 24) # 3
В целом, результаты не сильно изменились в сравнении с результатами без дамми-колонок, однако привели к изменению порядка состава главных компонент. В Dim-1 новая дамми-переменная ‘is Africa’ стала 6-й по вкладу в вариативность (около 7%), а новая дамми-переменная ‘is Oceania’ сильно не повлияла (во все три Dim ее влияние около 0%), что связано с большим числом (большей вариативности) объектов из Африки. Также, если вариативность дамми-переменных была бы значительна по сравнению со всеми остальными и было бы необходимо включить такие дамми-переменные в анализ, можно было бы стандартизировать все остальные переменные, чтобы дамми-переменные не искажали масштаб. В нашем случае целью было исследовать взаимосвязи между количественными переменными при помощи евклидовых расстояний, поэтому включение переменных ‘is Africa’ и ‘is Oceania’, которые по факту являются категориальными не лучшее решение.